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Abstract—A transient, numerical model for the prediction of the ground temperature at various depths 
below buildings is presented in this paper. The proposed model was developed by calculating the heat 
flow to the ground from a building, which depends on the complicated three-dimensional thermal process 
in the ground. The main difficulties in obtaining manageable solutions of the heat flow problem were: 
The three-dimensionality of the thermal process, the strong temporal variability of the outdoor temperature 
as well as the large number of parameters involved in describing the building foundation geometry as 
well as the thermal insulation. The techniques of superposition and numerical analysis were used to cope 
with these difficulties. The model was validated against experimental data and it was found that it could 
accurately predict the ground temperature under a building. 


1. INTRODUCTION 

The demand for more precise methods to calcu¬ 
late heating and cooling load of a building has 
been increased during the last two decades and 
has led to many studies of the transient heat 
flow from a building foundation to the ground. 

The use of direct or indirect earth coupling 
techniques for buildings requires knowledge of 
the ground temperature profile below the build¬ 
ing foundation. The direct earth contact involves 
partial or total placing of the building envelope 
in direct contact with the soil while the indirect 
contact involves the use of an earth-to-air heat 
exchanger system, (Mihalakakou et al., 1994a; 
Mihalakakou et al., 1994b; Santamouris et al, 
1995), through which air from the building or 
from the outside is circulated and then is 
brought into the building. 

The different calculation models, which were 
developed to predict the ground temperature 
under a building, are based on analytical solu¬ 
tions, on numerical analysis and experimental 
measurements. An important number of papers 
deal with a two-dimensional steady state heat 
flow process (Billington, 1951; Vuorelainen, 
1960a,b; Landman and Delsante, 1987. More¬ 
over, an extensive analysis of the heat flow 
processes in order to obtain sufficiently accurate, 
analytical and as simple as possible formulae, 
for the heat flow has been described by Claesson 
and Eftring (1980), Claesson and Dunand 
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(1983), Hagentoft (1988a), and by Claesson and 
Hagentoft (1991). 

Fourier analysis and transformation were 
used in Muncey and Spencer (1977) and in 
Delsante et al. (1983) to obtain expressions for 
the heat from uninsulated slabs. The boundary 
temperatures varied periodically with time while 
the temperature change under the walls, at the 
level of ground, was assumed to be linear. 

Delsante (1988) derived an analytical expres¬ 
sion for the two-dimensional steady state heat 
loss from a slab on the ground with the same 
thermal resistance at the ground surface and 
the floor. The temperature under the wall at the 
ground level was assumed to change linearly. 
A similar technique has been described in 
Hagentoft (1988b) and in Delsante (1989) while 
a method for calculating the heat loss of cellars 
was developed by Kusuda and Bean (1984). 

The temperature under a building foundation 
was numerically calculated for steady state cases 
by Adamson et al. (1973) and Backstrom et al. 
(1985). 

Boileau and Latta (1968) developed a simpli¬ 
fied method to calculate heat flow for uninsu¬ 
lated cellars. The results of this simplified 
method were compared with numerical predic¬ 
tions by Shipp and Broderick (1980) and it was 
found in a very good agreement for uninsulated 
building foundation. 

Mitalas (1982) developed a two-dimensional 
model for calculating the heat flow process for 
cellars. He divided the cellar wall and floor in 
five segments assuming that the outdoor temper¬ 
ature varied periodically. 
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Kusuda et al. (1983) described a calculation 
model in which the ground surface was divided 
in rectangular segments. The segment’s temper¬ 
ature varied periodically. The thermal problem 
was solved integrating the Green’s function 
solution. 

However, a large number of these methods 
are characterised by limited applicability as 
they do not take into account the three- 
dimensionality of the heat flow thermal process 
as well as the strong temporal variation of the 
outdoor temperature. 

The present paper deals with the rather com¬ 
plicated thermal process in the ground under a 
building and, in particular, with the heat flow 
through the foundation. The main goal was to 
obtain an accurate method to calculate the heat 
flow to the ground and thus the ground temper¬ 
ature below a building foundation at various 
depths. The present method can be used for a 
wide range of applications, such as the direct 
contact buildings as well as the indirect contact 
technique, and it can be easily used for the 
accurate prediction of the ground temperature 
profile when an earth-to-air heat exchangers’ 
system has been buried under a building. In this 
case, the ground temperature at a certain point 
in the pipe vicinity should be estimated by the 
superposition of the ground temperature caused 
by the heat exchanger presence and of the 
ground temperature caused by the ground 
surface boundary conditions. 

2. MODELLING OF GROUND 
TEMPERATURE UNDER A BUILDING 

The thermal problem of the heat flow from 
a building to the ground presents various 
difficulties caused by the genuinely three- 
dimensional, time-dependent character of 
the heat flow process in the ground below 
and around the building. The ground tempera¬ 
ture T(x,y,z,t) can accurately be modelled 
(Carslaw and Jaeger, 1959) by the following 
three-dimensional, transient, heat conduction 
equation: 

pC p (8T/8t) = \7(KV T) 

= d/dx(K8T/dx) + 8/8y{K8T/8y) 

+ 8/8z(K8T/8z) (1) 

where V indicates the gradient operator. 

The heat flux in the ground can be written as 
follows: 

-XV T= —K(8T/8x, 8T/8y , 8T/8z). (2) 


The total heat flow Q(t) from the building to 
the ground can be stated as follows: 

e,,) =JI -KVTndS. (3) 

Here n denotes the unit normal to S, pointing 
from ground to building. 

The initial condition is the following: 

T(x, y,z,t = 0) = T s (z) (4) 

where T s (z) was considered to be the undisturbed 
ground temperature. In order to solve the 
differential eqn (1) we assume that at t = 0 there 
is not any influence from the building’s presence. 
The boundary conditions are: 

(1) Co-ordinate x 
at x=x a : 

T(x=x„y,z,t)=T n {z) (5) 

at x=x b : 

T(x = x b ,y,z,t)=T u (z) (6) 

where, x a and x b are two large axial distances 
from the building foundation which were taken 
equal to 100 m. At those distances (x a , x b ) the 
temperature distribution is not influenced by 
the heat flow caused by the building presence. 
Thus, the temperature at x a and x b is the 
undisturbed soil temperature. T u (z) is the undis¬ 
turbed ground temperature at a certain depth z 
which does not change with time and is not 
influenced by the building’s presence. 

(2) Co-ordinate y 
at y=y a : 

T(x,y = y a ,z,t)=T u (z) (7) 

at y=y b : 

T(x, y=y h , z, t) = T u (z). (8) 

Similarly, y a and y b are two large distances in 
the y-direction where the temperature distribu¬ 
tion is not influenced by the building’s presence. 

(3) Co-ordinate z 

at z = 0 (at the ground surface): 

T(x,y,z = 0,r)=T boundary (9) 

where, T b ound ary was equal to the constant indoor 
surface temperature (T { ) for the case of the 
foundation interior region while for the case of 
the ground surface T boundary was equal to the 
outdoor temperature (T out (t)). 
at z = z c : 

T(x, y, z = z c , t) = T a (z c ) (10) 

where, z c was a large vertical distance in the 
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soil, taken equal to 50 m, where the soil temper¬ 
ature was considered to be the undisturbed soil 
temperature. 

Furthermore, the superposition technique was 
used to analyse the outdoor temperature varia¬ 
tion. Thus, the principle of superposition for 
this case was illustrated in Fig. 1 and for two 
superimposed components (1) and (2). 

The outdoor temperature at the ground sur¬ 
face can be represented (Carslaw and Jaeger, 
1959; Mihalakakou et al., 1992) by a Fourier 
series: 

T oul (t) = T 0 + | T n sin(2nnt/t p + 2. n ). (11) 

The two first harmonics of the eqn (11) have 
been used in this study in order to calculate the 
outdoor temperature at the ground surface. The 
thermal conductivity values used in this study 
were obtained from Puri (1986) and Gee (1966), 
while the soil density (p) and the soil specific 
heat capacity (C p ) were taken from Puri (1986). 
The transient building-ground system described 
in this paper incorporates four independent 
(x, y, z, t) and one dependent variables (T). The 
numerical method of control volume formula¬ 
tion was used to discretize the three-dimensional 
differential eqn (1). This method (Patankar, 
1980), which can be regarded as a special and 
new version of the method of weighted residuals, 
looks like the finite-difference method but it 
employs many ideas that are typical of the finite 
element methodology. 

The first step in the control volume formula¬ 


tion was to divide the calculation domain into 
a finite number of control volumes. A portion 
of a two-dimensional grid was shown in Fig. 2. 
For the grid point P, points E and W were its 
x-direction neighbours, while N and S are the 
y-direction neighbours. The control volume 
around P was shown by dashed lines. Finally, 
and for the present three-dimensional problem, 
two more neighbours T and B were added 
for the z-direction to complete the three- 
dimensional configuration. 

Moreover, it has been considered that the grid 
points were located at the geometric centre of 
each control volume. This location was pre¬ 
sented in Fig. 3 for a two-dimensional problem. 
The discretization equation was derived by 
integrating the differential eqn (1) over each 
control volume and over the time interval from 
t to r-F At. The order of this integration was 
chosen according to the nature of the term. For 
the representation of the terms AT/At was 
assumed that the grid point values of T prevail 
through the control volume. As it regards the 
terms AT/Ax, AT/Ay and AT/Az, it was neces¬ 
sary to make a profile assumption or an interpo¬ 
lation formula. The simplest possibility was to 
assume that the value of T at a grid point prevails 
over the control volume surrounding it. For this 
profile, the slope AT/Ax was not defined at the 
control volume faces. A profile which does not 
suffer from this difficulty is the piece-wise linear 
profile. In this study, linear interpolation func¬ 
tions were used between grid points. 

The time dependency was best handled 





Fig. 1. The technique of superposition where the temperature at the ground surface was divided in two 
components (1) and (2). 
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Fig. 2. Control volume for 

using implicit integration techniques. The fully 
implicit scheme satisfies the requirements of 
simplicity and physically satisfactory behaviour. 
Therefore, the discretization equation can easily 
be seen after the integration of the differential 
eqn (1) as follows: 

a P 7^, = a E T E + a w 7w + oc N T N + a s^s (12a) 


+ a x T t + a B T B + b 

where: 

a E = K e AyAz/(8x) e (12b) 

a w = K w AyAz/(dx) w (12c) 

a N = K n AzAx/(<5y) n (12d) 

a s = K s AzAx/(<5y) s (12e) 

a T = K t AxAy/(8z\ (12f) 

aB = K b AxAy/(Sz) b (12g) 

otp = pC p AxAyAz/At (12h) 

b = a$T> (12i) 


a P = a E + a w + a N + a s + a T + a B + ap. (12j) 


two-dimensional situation. 

The neighbours coefficients a E , oc w > oc N ,..., a B 
represent the conductance between the point P 
and the corresponding neighbour while the 
term a P Tp is the internal energy contained in 
the control volume at time t. 

In this study the whole system of the algebraic 
equations presented by eqn (12) were solved 
using the Gauss-Seidel iterative method, 
according to this method the values of the 
variables were calculated using each grid point 
in a certain order. 

This numerical method can be regarded as 
the best way to handle the serious difficulty of 
the large number of parameters involved in 
describing the foundation geometry and the 
thermal insulation. Therefore, considering a big 
number of control-volumes it was easier to 
specify the size and the shape of the foundation 
surface, the thickness of a thermal insulation 
and its thermal conductivity. 

The whole program was developed inside 
TRNSYS environment. TRNSYS (1990), is a 
transient system simulation program with a 
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O Control volume 


Fig. 3. Location of the control volume faces. 


modular structure which facilitates the addition 
to the program of other mathematical models 
not included in the standard TRNSYS library. 
The time step for the TRNSYS simulations was 
taken equal to 0.1 h. 

3. MODEL VALIDATION 

The collection and analysis of a comprehen¬ 
sive set of experimental data was considered an 
essential part of this research in order to provide 
a basis for the verification of the accuracy of 
the mathematical model. 

For this purpose, two experiments were car¬ 
ried out during the summer. Each experiment 
lasted 20 days, 11-30 June and 12-31 August. 
Data of the soil temperature were measured at a 
depth of 0.3 m below the earth surface and under 
the building foundation of the Philosophical 
Faculty of University of Ioannina. The floor 
area of this building was 2500 m 2 . The floor 
and the external walls were made of 60 mm 
dense concrete, 40 mm expanded polystyrene 
insulation and 60 mm dense concrete. The den¬ 
sity of dense concrete was equal to 2100 kg/m 3 


while the density of the polystyrene was equal 
to 25 kg/m 3 . The specific heat capacity of dense 
concrete was 833 J/kg°C while the specific 
heat capacity of the insulation was equal to 
1014J/kg°C. Accordingly, the thermal conduc¬ 
tivity of the dense concrete was 1.40 W/m°C 
and the conductivity of the insulation was 
0.03 W/m°C. The soil temperature data were 
recorded at 5 min intervals throughout the 
experiment. The ambient air temperature has 
been fluctuated during the experiments between 
9.8 and 39.3°C. 

The results of the experiments were compared 
with the theoretical calculations. Figures 4 and 
5 show the temporal variations of the measured 
and of the calculated soil temperature for the 
time periods 11-30 June and 12-31 August, 
respectively. As it can be seen from these figures 
there is a very good agreement between the 
predicted and measured values. Moreover, 
Fig. 6 compares the observed with the calculated 
soil temperature values for the whole time 
period of the two experiments. As shown, there 
is a very good agreement between the theoretical 
and the measured data while their maximum 
rarely exceeds 0.3°C. 
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Measured temperature values (°C) 

Fig. 6. Comparison between the observed with the calculated soil temperature values for the whole time period. 


CONCLUDING REMARKS 

An accurate, transient, implicit numerical 
model has been developed to calculate the 
ground temperature at various depths under 
a building foundation. The basic, three- 
dimensional thermal problem concerns homo¬ 
geneous ground, constant indoor temperature 
and an outdoor temperature presenting a large 
temporal variability. The proposed model was 
validated against experimental data and it was 
found that it could accurately predict the 
temperature of the soil under a building 
foundation. 


NOMENCLATURE 

C p specific heat capacity (J/kg°C) 

K soil thermal conductivity (W/m°C) 

K, the thermal conductivity corresponding to grid 
point I (W/m°C) 

Q heat flow to the ground (W) 

Sj notation for the surface foundation against the 
ground (m 2 ) 
t time (s) 

f p time period for outdoor periodic temperature (s) 

T temperature in the soil (°C) 

^boundary temperature at the boundary (°C) 

7j indoor temperature (°C) 

T l the temperature of the grid L (°C) 

T n amplitude of the temperature wave at the ground 
surface (°C) 


T 0 annual mean temperature at the ground surface 
(°C) 

T out (f) outdoor temperature (°C) 

T p ° the old value (at time t) of the variable (°C) 

T u (z) the undisturbed soil temperature, which is not 
influenced by the building presence (°C) 
x,y horizontal Cartesian co-ordinates (m) 
z vertical Cartesian co-ordinate (m) 

Greek characters 

Ax x-direction width of the control volume (m) 

Ay y-direction width of the control volume (m) 

Ar z-direction width of the control volume (m) 

Sx x-direction distance between two adjacent grid 
points (m) 

5y y-direction distance between two adjacent grid 
points (m) 

Sz z-direction distance between two adjacent grid 
points (m) 

p soil density (kg/m 3 ) 

A n phase for periodic temperature 

Subscripts 

B neighbour in the z-direction (bottom) 
b control-volume face between P and B 
E neighbour in the positive x-direction (east side) 
e control-volume face between P and E 
N neighbour in the positive y-direction (north side) 
n control-volume face between P and N 
P central grid point under consideration 
S neighbour in the negative y-direction (south side) 
s control-volume face between P and S 
T neighbour in the positive z-direction (top) 
t control-volume face between P and T 
W neighbour in the negative x-direction (west side) 
w control-volume face between P and W 
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